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Abstract 

The corrections to finite-size scaling in the critical two-point correlation 
function G(r) of 2D Ising model on a square lattice have been studied numer- 
ically by means of exact transfer-matrix algorithms. The systems have been 
considered, including up to 800 spins. The calculation of G(r) at a distance 
r equal to the half of the system size L shows the existence of an amplitude 
correction cx L~ 2 . A nontrivial correction cx L~ - 25 of a very small magni- 
tude also has been detected, as it can be expected from our recently developed 
GFD (grouping of Feynman diagrams) theory. Monte Carlo simulations of the 
squared magnetization of 3D Ising model have been performed by Wolff's algo- 
rithm in the range of the reduced temperatures t > 0.000086 and system sizes 
L < 410. The effective critical exponent (3 e ff{t) tends to increase above the 
currently accepted numerical values. The critical coupling K c — 0.22165386(51) 
has been extracted from the Binder cumulant data within L 6 [96; 384]. The 
critical exponent l/v, estimated from the finite-size scaling of the derivatives 
of the Binder cumulant, tends to decrease slightly below the RG value 1.587 for 
the largest system sizes. The finite-size scaling of accurately simulated maximal 
values of the specific heat Cy in 3D Ising model confirms a logarithmic rather 
than power-like critical singularity of Cy. 

Keywords: Transfer matrix, Ising model, if 4 model, critical exponents, finite- 
size scaling, Monte Carlo simulation 
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1 Introduction 

Since the exact solution of two-dimensional Lenz-Ising (or Ising) model has been 
found by Onsager fl_, a study of various phase transition models is of permanent 
interest. Nowadays, phase transitions and critical phenomena is one of the most 
widely investigated fields of physics [21 El- Remarkable progress has been reached 
in exact solution of two-dimensional models Hf. Recently, we have proposed [5| a 
novel method based on grouping of Feynman diagrams (GFD) in ip 4 model. Our 
GFD theory allows to analyze the asymptotic solution for the two-point correlation 
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function at and near criticality, not cutting the perturbation series. As a result the 
possible values of exact critical exponents have been proposed jS] for the Ginzburg- 
Landau (p 4 ) model with 0(n) symmetry, where n = 1, 2, 3, ... is the dimensionality 
of the order parameter. Our predictions completely agree with the known exact 
and rigorous results in two dimensions [I]. In [3], we have compared our results 
to some Monte Carlo (MC) simulations and experiments [HI El El- The examples 
considered there support our predictions about the critical exponents. A more recent 
comparison with experimental data very close to the A-transition point T = T\ in 
liquid helium has been made in |2j . We have shown there that our critical exponents 
better describe the closest to T\ data for the superfluid fraction of liquid helium as 
compared to the exponents provided by the perturbative renormalization group (RG) 
theory |1U1 1111 112j . As claimed in [3] , the conventional RG expansions are not valid 
from the mathematical point of view. The current paper, dealing with numerical 
transfer-matrix analysis of the two-point correlation function in 2D Ising model, 
as well as with MC simulations in the three-dimensional Ising model presents some 
more evidences in favour of the critical exponents predicted by the GFD theory. Our 
estimations are based on the finite-size scaling theory, which by itself is an attractive 
field of investigations ^Hj and has increasing importance in modern physics [H]. 



2 Critical exponents predicted by GFD theory 

Our theory predicts possible values of exact critical exponents 7 and v for the ip 4 
model whith 0(n) symmetry (n-component vector model) given by the Hamiltonian 

H/T = J [r p\x) + c(V<^(x)) 2 + V(x)] dx , (1) 

where tq is the only parameter depending on temperature T, and the dependence 
is linear. At the spatial dimensionality d = 2, 3 and n = 1,2,3,... the predicted 
possible values of the critical exponents are [Hj 

d + 2j + 4m 

7 d{l + m + j) - 2j ' U 



2(1 + m) + j 
d(l + m + j)- 2j 



(3) 



where m > 1 and j > —m are integers. It is well known that the 0(n)-symmetric 
ip 4 model belongs to the same universality class as the corresponding lattice model 
(Ising model at n = 1, XY model at n = 2, the classical Heisenberg model at 
n = 3, etc.), where the order parameter is an n-component vector (spin) with fixed 
modulus j </?(x) |= 1, since the latter is a particular case (ro — ► — 00 at — ro/(2u) = 1 
or A — > 00 in the notations used in ^lj) of the lattice ip 4 model, where the gradient 
term is represented by finite differences [14] . Besides, the partition functions and 
two-point correlation functions of both (p 4 model in [5] and Ising model can be 
represented by similar functional integrals [15] 116] . Thus, at n = 1 we have m = 3 
and j = to fit the known exact results for the two-dimensional Ising model. As 
proposed in Ref. j^j, in the case of n = 2 we have m = 3 and j = 1, which yields in 
three dimensions v = 9/13 and 7 = 17/13. 

As already explained in [3], our predictions do not refer to the case of the self- 
avoiding random walk recovered at n = 0. The values (J2J) and (jSJ) have been derived 
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in [3] assuming that 2^ — 7 > holds. In principle, the mean-field-like solution 
with 2v — 7 = can exist at d < 4, and it refers to the Gaussian random walk with 
n = —2. This is a special case, not related to (J2J and (JHJ, where two expansion 
parameters A 2l/_7 and A 27_cft/ with A = T — T c — > being the deviation from the 
critical temperature T c , are replaced by one parameter A 27_dl/ . Eq. (48) in [5] can 
be then satisfied with 7 = 1 and v = 1/2, i. e., all the exponents are consistent 
and each term can be compensated. The singularity of the specific heat with the 
exponent a = 2 — d/2 comes from the leading terms in Eq. (60) of jSJ. Obviously, 
7 = 1 and v = 1/2 always are the true exponents at d > 4, where the Gaussian 
approximation G(k) = 1/ [G(0) _1 + 2ck 2 ] for the two-point correlation function 
G(k) in the Fourier representation is asymptotically exact at u — > and T > T c 
for arbitrarily small wave vectors k. These exponents are recovered at any m and j 
in © and |j3J) when approaching the upper critical dimension d = 4 from below. 

Our formulae do not provide any sensible result approaching d = 1, where v = 
l/(d — 1) is expected at n = 1 according to the MidgaPs approximation ^7]. It can 
be understood from the point of view that 2, probably, is the marginal value of 
d, such that an analytic continuation of the results from <i-dimensional hypercubes 
can be only formal and has no physical meaning at d < 2. In this sense, we expect 
that d = 2 is a special dimension for any n > 1. Besides, the critical temperature 
does not vanish at d — > 2 + 0, and for n = 1 there exist lattices for which the critical 
temperature is nonzero at the fractal dimension below 2 [T^]. In the marginal case 
d = 2 different behavior is observed at low temperatures: the long-range order at 
n=l, the Kosterlitz-Thouless structural order at n = 2, and disordered state at 
n > 2. 

Our concept agrees with the known rigorous results for XY model |2()1 I21j . 
It disagrees with the prediction of the perturbative RG theory [22] that the critical 
temperature goes to zero at d — > 2+0 for the 0{n > 2)-symmetric nonlinear a model 
and, therefore, the behavior in this case is Gaussian, i. e., rj = and v = l/(d — 2). 
The results of the perturbative RG theory are not rigorous since the claims are based 
on formal expansions which break down in relevant limits, in this case at vanishing 
external field H — > +0. Moreover, essential claims of this theory are based on an 
evidently incorrect mathematical treatment. In particular, the conclusion about 
the Gaussian character of the 0(n)-symmetric ip 4 model below T c has been made 
in by simply rewritting the Hamiltonian in an apparently Gaussian form (see 
Eqs. (3.4) to (3.6) in [2'6\). The author, however, forgot to include the determinant of 
the transformation Jacobian in the relevant functional integrals, according to which 
the resulting model all the same is not Gaussian. Due to the reasons mentioned 
above, we do not believe in predictions of the perturbative RG theory, but rely only 
on exact and rigorous results. 

There exists a simple non-perturbative explanation why the critical temperature 
should stay finite at d — ► 2 + for the 0(n)-symmetric Heisenberg model. Below T c , 
the difference in free energies for models with antiperiodic and periodic boundary 
conditions along one axis is AF oc T(T)L d ~ 2 , where L is the linear size of the 
system and T(T) is the helicity modulus. It holds because the energy difference in 
the ground state at T = is oc L d ~~ 2 , corresponding to gradually rotated spins in 
any given plane. The factor T(T) takes into account the temperature dependence. 
It vanishes at T > T c . Hence, the factor L d ~ 2 always vanishes at d < 2 in the 
thermodynamic limit L — > 00, therefore the long-range order (if it would exist) 
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could be destroyed in this case at any finite temperature by gradually rotating the 
spins without increasing the free energy. Thus, the long-range order disappears at 
d < 2 irrespective to the behavior of T(T), i. e., irrespective to the value of T c 
at d = 2 + 0. In such a way, the assumption that the critical temperature should 
go to zero continuously appears as an additional unnecessary constraint. On the 
other hand, if the critical temperature remains finite in (£> 4 model, then rj should be 
positive at d = 2 + to avoid the divergence of (ip 2 (x.)} = n L~ d J2k G(\z) at T = T c , 
where G(k) ~ a k~ 2+v with a / holds for the two-point correlation function. The 
expectation rj > agrees with (J2J and ©• However, the critical temperature at 
d = 2 + and rj tend to zero in the limit n — > oo to coincide with the known exact 
results for the spherical model, which are recovered in ((2j) and © at j/m — > oo. 

In the present analysis the correction-to-scaling exponent 6 for the susceptibility 
is also relevant. The susceptibility is related to the correlation function G(k) via 
X oc G(0) [llj . In the thermodynamic limit, this relation makes sense at T > T c . 
According to our theory, i 7 G(0) can be expanded in a Taylor series of t 2u ~ 1 at t — > 0. 
In this case the reduced temperature t is defined as t = ro(T) — ro(T c ) oc T — T c . 
Formally, t 2 ~ l ~ dv appears as second expansion parameter in the derivations in Ref. [3], 
but, according to the final result represented by Eqs. (J2J) and ©, (27 — dv) / (2v — 7) 
is a natural number. Some of the expansion coefficients can be zero, so that in 
general we have 

e = e(2u- 7 ), (4) 

where t may have integer values 1, 2, 3, etc. One can expect that £ = 4 holds at 
n = 1 (which yields 9 = 1 at d = 2 and 9 = 1/3 at d = 3) and the only nonva- 
nishing corrections are those of the order t e , t 29 , t 39 , since the known corrections 
to scaling for physical quantities, such as magnetization or correlation length, are 
analytical in the case of the two-dimensional Ising model. Here we suppose that 
the confluent corrections become analytical, i. e. 9 takes the value 1, at d = 2. 
Besides, similar corrections to scaling are expected for susceptibility \ an d mag- 
netization M since both these quantities are related to G(0), i. e., x « G(0) and 
M 2 = \im x -, oo (ip{0)ip(x)) = limy^oo G(0)/V hold where V = L d is the volume and 
L is the linear size of the system. The above limit is meaningful at L — > 00, but 
G(0)/V may be considered as a definition of M 2 for finite systems too. The latter 
means that corrections to finite-size scaling for \ and M are similar at T = T c . Ac- 
cording to the scaling hypothesis and finite-size scaling theory, the same is true for 
the discussed here corrections at t — > 0, where in both cases (x and M) the defini- 
tion t =| J"o(T) — ro(T c ) I is valid. Thus, the expected expansion of the susceptibility 
X looks like x = t ~ 7 ( a o + a it 9 + Q^ 26 * + ■ • -J . In this discussion we have omitted 
the irrelevant for critical behaviour background term in the susceptibility, which is 
constant in the first approximation and comes from the short-distance contribution 
to x = J2*.G(x) JJ, where G(x) is the real-space two-point correlation function. 

Our hypothesis is that j = j(n) and I = £(n) monotoneously increase with n to 
fit the known exponents for the spherical model at n — > 00. The analysis of the MC 
and experimental results here and in [5] enables us to propose that j(n) = n — 1, 
£(n) = n + 3, and m = 3 hold at least at n = 1,2,3. These relations, probably, 
are true also at n > 3. This general hypothesis is consistent with the idea that the 
critical exponents 7, u, and 9 can be represented by some analytical functions of n 
which are valid for all natural positive n and yield rj = 2 — 7/z^ oc 1/n rather than 
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i] oc l/n s with s = 2,3, . . . (s must be a natural number to avoid a contradiction, 
i. e., irrational values of j{n) at natural n) at n — ► oo. At these conditions, j( n ) 
and are linear functions of n (with integer coefficients) such that £(n)/j(n) — > f 
at n — » oo, and m is constant. Besides, j(l) = 0, m(l) = 3, and ^(1) = 4 hold 
to coincide with the known results at n = 1. Then, our specific choice is the best 
one among few possibilities providing more or less reasonable agreement with the 
actually discussed numerical an experimental results. 

We allow that different i values correspond to the leading correction-to-scaling 
exponent for different quantities related to G(k). According to [jj], the expansion of 
Cr(k) in <f model by itself contains a nonvanishing term of order i 2i/_7 = (in the 
form G(k) ~ t~T ^(kt - ") + t^ v g x (kt~ v )] whith 5l (0) = 0, since £ > 1 holds in the 
case of susceptibility) to compensate the corresponding correction term (produced 
by c(V</?) 2 ) in the equation for 1/G(k) (cf. |S]). 

The correction is related to correction term in the long-distance (x — ► oo) 
behavior of the real-space pair correlation function G(x) oc x 2 ~ d ~ v [1 + (x~ v )] at 
the critical point, as well as to correction in the finite-size scaling expressions 
at criticality. Such kind of corrections must not necessarily appear in the Ising 
model, where they could have zero amplitude. In particular, the critical Green's 
(correlation) function G(x) of 2D Ising model in (11) crystallographic direction on 
an infinite lattice can be calculated easily based on the known exact formulae 24 , 
and it yields G(x) oc x -1 / 4 [1 + (^~ 2 )] at large distances x — * oo. Nevertheless, 
our calculations in 2D Ising model discussed in Sec. 14.31 indicate the existence of 
a nontrivial finite-size correction of the kind L~ v (for (10) direction), as it can be 
expected from our theoretical results for the ip 4 model. The thermodynamic limit 
is a particular case of the finite-size scaling with the scaling argument x/L — * 0, 
therefore it is possible that the nontrivial corrections to the correlation function in 
2D Ising model vanish in this special case or in the crystallographic direction (11), 
but not in general. 

Our consideration can be generalized easily to the case where the Hamiltonian 
parameter tq is a nonlinear analytical function of T. Nothing is changed in the 
above expansions if the reduced temperature t, as before, is defined by t = ro(T) — 
ro(T c ). However, analytical corrections to scaling appear (and also corrections like 
(T - T c ) m+nd with inte ger m and n) if t is reexpanded in terms of T — T c at T > T c . 
The solution at the critical point remains unchanged, since the phase transition 
occurs at the same (critical) value of r$. 

3 Exact transfer matrix algorithms for calculation of 
the correlation function in 2D Ising model 

3.1 Adoption of standard methods 

The transfer matrix method, applied to analytical calculations on two-dimensional 
lattices, is well known PJI1]. The asymptotic behavior of the correlation functions 
can be studied by means of the equations of the conformal field theory Exact 
equations for the two-point correlation function of 2D Ising model on an infinite 
lattice are known, too )24| . However, no analytical methods exist for an exact 
calculation of the correlation function in 2D Ising model on finite-size lattices. This 
can be done numerically by adopting the conventional transfer matrix method and 
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Figure 1: Illustrative examples of the lattices with dimensions NxL (a) and V2N x V2L (b) 
with periodic boundary conditions along the dashed lines. The correlation function has been 
calculated in the (10) crystallographic direction, as indicated by the arrows. 



modifying it to reach the maximal result (calculation of as far as possible larger 
system) with minimal number of arithmetic operations, as discussed further on. 

We consider the two-dimensional Ising model where spins are located either on 
the lattice of dimensions NxL, illustrated in Fig. or on the lattice of dimensions 
y/2N x \/2L, shown in Fig. ^p. The periodic boundaries are indicated by dashed 
lines. In case (a) we have L rows, and in case (b) - 2L rows, each containing N 
spins. Fig. Q shows an illustrative example with N = 4 and L = 3. In our notation 
nodes are numbered sequently from left to right, and rows - from bottom to top. 

For convenience, first we consider an application of the transfer matrix method 
to calculation of the partition function 

Z=£exp Wj , (5) 

Wk} V (id) I 

where <7j = ±1 are the spin variables, and the summation runs over all the possible 
spin configurations {crfc}. The argument of the exponent represents the Hamiltonian 
of the system including summation over all the neighbouring spin pairs of the 
given configuration {cr^}; parameter (3 is the coupling constant. Let us consider 
lattice (a) in Fig. ^ but containing n rows without periodic boundaries along the 
vertical axis and without interaction between spins in the upper row. We define 
the 2 Ar -component vector r n such that the i— th component of this vector represents 
the contribution to the partition function provided by the i— th spin configuration 
of the upper row. Then we have a recurrence relation r n+ i = Tr n , where T is 
the transfer matrix which includes the Boltzmann weights of newly added bonds. 

Furthermore, we can write r^_-^ — T^r^, where is the partial contribution to r n 

(i) 

provided by the i-th configuration of the first row. The components of are given 
by j = <5j,i- In the case of the periodic boundary conditions the (L + l)-th row 
must be identical to the first one, which leads to the well known expression [H I26j 

^ = E( r ai) 4 = ^-e( rL )=E^ (6) 

i i 
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where Aj are the eigenvalues of the transfer matrix T. An analogous expression for 
the lattice in Fig. reads 

^ = E( r S + i) j = Trace([T 2 T 1 ] L ) , (7) 

i 

where the vectors obey the reccurence relation = Vn with different 

transfer matrices T\ and T2 for odd and even row numbers n, respectively. 

The actual scheme can be easily adopted to calculate the correlation function 
(fTjCTj) between any two lattice points i and j. Namely, the correlation funtion G(x) 
between the points separated by a distance x, like indicated in Fig. Q is given by 
the statistical average Z'/Z, where the sum Z' is calculated in the same way as Z, 
but including the corresponding product of spin variables. It implies the following 
replacements: 

(r?) . = 5 hi =* (if) . = 5 hl (n- 1 jr [a^ [a(£ + x)]\ : case (a) (8) 

3 3 V i=i J 

(4) . =* (rfli) . x [ N- 1 J2 W% [a(£ + A(x))]\ : case (b) , (9) 

where A(x) = x/2 holds for even x, and A(x) = (x — l)/2 — for odd x. In our 
notation, [cr(/c)]j is the spin variable in the fc-th node in a row provided that the 
whole set of spin variables of this row forms the i-th configuration. It is supposed 
that a(k + N) = a(k) holds according to the periodic boundary conditions. Such a 
symmetrical form, which includes an averaging over I, allows to reduce the amount 
of numerical calculations: due to the symmetry we need the summation over only 
~ 2 N /N nonequivalent configurations of the first row instead of the total number of 
2^ configurations. 



3.2 Improved algorithms 

The number of the required arithmetic operations can be further reduced if the 
recurrence relations r^ +1 = Tr„' and r^ +l = Ti^rff are split into N steps of 
adding single spin. To formulate this in a suitable way, let us first number all the 
2^ spin configurations {<r(l); c(2); • • • ; a(N — 1); cr(N)} by an index i as follows: 



i 


= 1 : 


{-i; 


-1; ••• 


-i; 


-i; 


-1} 


i 


= 2 : 


{-i; 


-1; ••• 


-i; 


-i; 


+1} 


i 


= 3 : 


{-i; 


-1; ••• 


-i; 


+i; 


-1} 


i 


= 4 : 


{-i; 


-1; ••• 


-i; 


+i; 


+1} 


i 


= 2 N : 




+1; ••• 


+i; 


+i; 


+1} 



(10) 



We remind that the sequence of numbers in the i-th row corresponds to the spin 
variables [cr(/c)]j with k = 1,2, . . . , N . They change with i just like the digits of 
subsequent integer numbers in the binary counting system. 

Consider now a lattice where n rows are completed, while the (n + l)-th row 
contains only I spins where I < N, as illustrated in Fig. [21 in both cases (a) and (b) 
taking as an example N = 4. We consider the partial contribution (r n+ i (i. e., 
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Figure 2: Schematic pictures illustrating the algorithms of calculation for the lattices a and 
b introduced in Fig. ^ 

i—th component of vector r n+ i^) in the partition function Z (or Z') provided by 
a fixed (i—th) configuration of the set of N upper spins. These are the sequently 
numbered spins shown in Fig.[21by empty circles. For simplicity, we have droped the 
index denoting the configuration of the first row. In case (b), the spin depicted by 
a double-circle has a fixed value a'. In general, this spin is the nearest bottom-left 
neighbour of the first spin in the upper row. According to this, one has to distinguish 
between odd and even n: a' refers either to the first (for odd n), or to the iV-th 
(for even n) spin of the n-th row. It is supposed that the Boltzmann weights are 
included corresponding to the solid lines in Fig. [2] connecting the spins. In case (a) 
the weights responsible for the interaction between the upper numbered spins are 
not included. Obviously, for a given £ > 1, r n+ i t £ can be calculeted from r n +i^_i via 
summation over one spin variable, marked in Fig. Elby a cross. In case (a) it is true 
also for £ = 1, whereas in case (b) this variable has fixed value a 1 at I = 1. In the 
latter case the summation over a' is performed at the last step when the (n + l)-th 
row is already completed. These manipulations enable us to write 

r n+1 = Tr n = W N W N -i ■■■W 2 W 1 r n (11) 

with 

w t =Y, > ( 12 ) 

a=±l 

where the componets of the matrices We(cr) are given by 

(WiCa)).. = 5(j,ii(^l^))-exp(/3a{[ C T(l)] J + [a(2)] l + [ CT (iV)]J) 
(W £ (a)) tJ = SVJiiaAty'exptfviWtyi + We+lM) : 1 <£< N 
(W N (a)) i:j = 5(i,i 1 ((T ) iV,i))-exp(/3 ( 7[ ( r(iV)] i ) . (13) 

Here 5(j, k) is the Kronecker symbol and 

j 1 (a,l,i)=i + (a-[a(£)] l )2 N - e ' 1 (14) 

are the indexes of the old configurations containing £—1 spins in the (n + 1)— th row 
depending on the value a of the spin marked in Fig. by a cross, as well as on the 
index i of the new configuration with £ spins in the (n + l)-th row, as consistent 
with the numbering (|1()[). 
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The above equations (fTT|) to (|T3|) refers to ). In case (b) we have 

r n+1 = T 1>2 r n = W> W§» ■ ■ ■ (a') r n , (15) 

cr'=±l 

— (1 2) 

where W) ' are the matrices 

W? A = E ^fV) . (16) 

cr=±l 

Here indexes 1 and 2 refer to odd and even row numbers n, respectively, and the 
components of the matrices W) ' (a) are 

(W^ 2 \aj) .. = 5 (jJiflM 0) • exp (/? [a(£)] t {a + [a(t + 1)].}) , (17) 

where [<r(iV + l)]j = <r' and the index ji(<r, i) is given by (|14|). For the other index 
we have 

J2 (a,l,i) = 2i-2^ 1 ([ f r(l)] t + l) + i(a-l) 

i 2 («r,A») = : ^>2. (18) 

— ' — '(1 2) 

Note that the matrices Wg and ' have only two nonzero elements in each 
row, so that the number of the arithmetic operations required for the construction 
of one row of spins via subsequent calculation of the vectors r n+ i ^ increases like 
2N ■ 2 N instead of 2 2N operations necassary for a straightforward calculation of the 
vector Tr n . Taking into account the above discussed symmetry of the first row, the 
computation time is proportional to 2 2L L for both L x L (a) and y/2L x y/2L (b) 
lattices in Fig. ^ with periodic boundary conditions. 

3.3 Application to different boundary conditions 

The developed algorithms can be easily extended to the lattices with antiperiodic 
boundary conditions. The latter implies that a(N + 1) = —a(N) holds for each 
row, and similar condition is true for each column. We can consider also the mixed 
boundary conditions: periodic along the horizontal axis and antiperiodic along the 
vertical one, or vice versa. To replace the periodic boundary conditions with the 
antiperiodic ones we need only to change the sign of the corresponding products of 
the spin variables on the boundaries. Consider, e. g., the case (a) in Fig. ^ The 
change of the boundary conditions along the vertical axis means that the first term 
in the argument of the exponent in each of the Eqs. (|13|) changes the sign for the last 
row, i. e., when n = L. The same along the horizontal axis implies that the term 
[iy(N)] i in the equation for {Wi{a))- changes the sign. In this case, however, the 
symmetry with respect to the configurations of the first row is partly broken and, 
therefore, we need summation over a larger number of nonequivalent configurations. 

4 Transfer matrix study of critical Greens function and 
corrections to scaling in 2D Ising model 

4.1 General scaling arguments 

It is well known that in the thermodynamic limit the real-space Greens function of 
the Ising model behaves like G(r) oc r 2 ~ d ~ v at large distances r — ► oo at the critical 
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point [5 = /3 C , where r\ is the critical exponent having the value rj = 1/4 in two 
dimensions (d = 2). Based on our transfer matrix algorithms developed in Sec 01 
here we test the finite-size scaling and, particularly, the corrections to scaling in 2D 
Ising model at the critical point (3 = (3 C = \ In (l + V2). 

In [S] the critical correlation function in the Fourier representation, i. e. G(k) at 
T = T c , has been considered for the model. In this case the minimal value of the 
wave vector magnitude k is related to the linear system size L via k m i n = 2tt/L. In 
analogy to the consideration in Sec. 5.2 of [3], one expects that k/k m i n is an essential 
finite-size scaling argument, corresponding to r/L in the real space. In the Ising 
model at r ~ L one has to take into account also the anisotropy effects, so that the 
expected finite-size scaling relation for the real-space Greens function at the critical 
point (3 = f3 c reads 

G(r) ~r 2 - d -' ri f(r/L) : r^oo,L^oo, (19) 

where the scaling function f{z) depends also on the crystallographic orientation of 
the line connecting the correlating spins, as well as on the orientation of the periodic 
boundaries. A natural extension of (|19|) . including the corrections to scaling, is 

G(r) = J2r- X *fi(r/L), (20) 

where the term with Ao = d — 2 + r/ is the leading one, whereas those with the 
subsequently increasing exponents Ai, A2, etc., represent the corrections to scaling. 
By a substitution fe(z) = z Xe f'(z), the asymptotic expansion (|20j) transforms to 

G(r) = f>(r/L) L~ x ° (l + £ L~»' f t (r/L) j , (21) 

where fe(z) = fi(z)/fo(z) and tug = Xp — Ao are the correction-to-scaling exponents. 

We have tested the scaling relation (|19|) in 2D Ising model by using the exact 
transfer matrix algorithms in Sec. 01 We have found that all points of f{r/L) = 
r 1 / 4 G(r) for the correlation function in (10) direction (case (a) in Sec. OJ well fit a 
common smooth line at 2 < r < L/2 and L = 8, 12, 15, and 18. It implies that the 
corrections to ()19|) are rather small. 

4.2 Correction— to— scaling analysis for the L x L lattice 

Based on the scaling analysis in Sec. 14.11 here we discuss the corrections to scaling 
for the lattice in Fig. We have calculated the correlation function G(r) at a fixed 
ratio r/L = 0.5 in (10) direction, as well as at r/L = 0.5\/2 in (11) direction at 
L = 2, 4, 6, . . . with an aim to identify the correction exponents in 1)21(1 . Note that in 
the latter case the replacement @ is valid for G{\/2x) (where x = 1, 2, 3, . . .) with 
the only difference that A(x) = x. 

Let us define the effective correction-to-scaling exponent uj e ff(L) in 2D Ising 
model via the solution of the equations 

L x / 4 G(r = const -L) = a + b L~"*ff (22) 
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Table 1: The critical correlation function G(r — c ■ L) in (10) (c = 0.5) and (11) (c = 
0.5V2) crystallographic directions vs the linear size L of the lattice (a) in Fig. and the 
corresponding effective exponents u e ff(L) and uj(L). 



L 


direction (10) 


direction (11) 


G(0.5L) 


Ueff(L) 


G(0.5\/2L) 


We// (£) 


u(L) 


2 


0.84852813742386 


2.7366493 


0.8 


1.8672201 




4 


0.74052044609665 


2.9569864 


0.71375464684015 


2.2148707 




6 


0.67202206468538 


1.8998036 


0.65238484475089 


2.1252078 




8 


0.62605120856389 


1.5758895 


0.60935351016910 


2.0611362 


1.909677 


10 


0.59238112628953 


1.6617494 


0.57724041054810 


2.0351831 


1.996735 


12 


0.56615525751968 


1.7774398 


0.55200680271678 


2.0232909 


2.002356 


14 


0.54485584658226 


1.8542943 


0.53141907668442 


2.0167606 


2.001630 


16 


0.52703456475995 




0.51414720882560 






18 


0.51178753041103 




0.49934511003360 







at L = L, L + AL, L + 2AL with respect to three unknown quantities u e ff, a, and 
b. According to (|2T|) . where Ao = rj = 1/4, such a definition gives us the leading 
correction-to-scaling exponent u> at L — > 00, i. e., lirn£_ ) . 00 uj e ff{L) = a;. 

The calculated values of G(r = c - L) in the (10) and (11) crystallographic direc- 
tions [in case (a)] with c = 0.5 and c = 0.5\/2, respectively, and the corresponding 
effective exponents u e ff(L), determined at AL = 2, are given in Tab. ^ I n both 
cases the effective exponent u e ff(L) seems to converge to a value about 2. Besides, 
in the second case the behavior is smoother, so that we can try someway to extrap- 
olate the obtained sequence of w e ff values (column 5 in Tab.^) to L = 00. For this 
purpose we have considered the ratio of two subsequent increments in uj e ff, 

r(T) = ^eff{L + AL)-UJ eff {L) 

K > u eff {L)-u efi {L-AL) ■ K ' 

A simple analysis shows that r(L) behaves like 

r(L) = 1 - AL ■ {(J + l)/,- 1 + O (l- 2 ) (24) 



at L — > 00 if uj e ff(L) =uj + holds with an exponent u>' > 1. The numerical 

data in Tab. ^ show that Eq. (|24[) represents a good approximation for the largest 
values of L at uj' = 2. It suggests us that the leading and the subleading correction 
exponents in (|2"T|) could be oj = u\ = 2 and u>2 = 4, respectively. Note that u e ff(L) 
can be defined with a shift in the argument. Our specific choice ensures the best 
approximation by Q24[) at the actual finite L values. 

Let us now assume that the values of u e f j (L) are known up to L = L max . Then 
we can calculate from (|23|) the r(L) values up to L = L max — AL and make a suitable 
ansatz like 

r(L) = 1 - 3AL- ZT 1 + 6ZT 2 at L > L max (25) 
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Table 2: The critical correlation function G(r = L) in (10) crystallographic direction and 
the effective exponents u) e tf(L) and u>{L) vs the linear size L of the lattice (b) in Fig.^ 



L 


G(L) 


Weff(L) 




2 
3 





0.8 

7203484812087670 










4 





6690636562097066 










5 





6321925914229602 










6 





6037455936471098 










7 


u 


580 /booo049zb8o8 










8 





5616046762441826 


2 


066235298 






9 





5452468033693456 


2 


043461090 






10 





5310294874153481 


2 


030235674 


1 


996772124 


11 





5184950262041604 


2 


022130104 


1 


999333324 


12 





5073151480587211 


2 


016864947 


1 


999941357 


13 





4972468711401118 


2 


013265826 


2 


000036957 


14 





4881056192765374 


2 


010701166 


2 


000040498 


15 





4797481011874659 


2 


008811505 


2 


000044005 


16 





4720609977942179 


2 


007380630 


2 


000053415 


17 





4649532511721054 


2 


006272191 


2 


000063984 


18 





4583506666254706 


2 


005396785 


2 


000073711 


19 





4521920457268738 










20 





4464263594840965 











for a formal extrapolation of uj e ff(L) to L = oo. This is consistent with (|24jl where 
a/ = 2. The coefficient b is found by matching the result to the precisely calculated 
value at L = L max — AL. The subsequent values of cu e //(L), calculated from (|2"H|) 
and $ZE§ at L > L max , converge to some value uj(L max ) at L — > oo. If the leading 
correction-to-scaling exponent w is 2, indeed, then the extrapolation result Lu(L max ) 
will tend to 2 at L max — * oo irrespective to the precise value of a/. 

As we see from Tab. [l\ the values of u;(L) come remarkably closer to 2 as com- 
pared to uj e ff(L), suggesting that u = 2. As we have discussed in Sec.|2 there could 
be a nontrivial correction in 1)21(1 with u = r] = 1/4. If it really exists, then it has a 
very small amplitude, otherwise it would show up in our analysis. 

4.3 Correction-to-scaling analysis for the \plL x ^/2\L lattice 

To test the possible existence of nontrivial corrections to scaling, here we make the 
analysis of the correlation function G(r) in (10) direction on the \f2L x y/2L lattice 
shown in Fig. ^>. The advantage of case (b) in Fig. ^ as compared to case (a) is 
that y/2 times larger lattice corresponds to the same number of the spins in one 
row. Besides, in this case we can use not only even, but all lattice sizes to evaluate 
the exponent uj from calculations of G(r = L), which means that it is reasonable 
to use the step AL = 1 to evaluate uJ e ff an d from Eqs. (|22[). and 1)25(1 . 

The results, are given in Tab. [2j It is evident from Tab. that the extrapolated 
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Figure 3: The deviation of the extrapolated effective exponent Auj(L) = ui(L) — 2 as a 
function of L~ 4 . The extrapolation has been made by using the calculated G(r) values in 
Tab. up to the size L + 2. A linear convergence to zero would be expected in absence of 
any correction term with exponent lu < 2. 

values of the effective correction exponent, i. e. ui(L), come surprisingly close to 2 
at certain L values. Besides, the ratio of increments r [cf. Eq. (|23|) ] in this case is 
well approximated by (|25[). as consistent with existence of a correction term in (|21l) 
with exponent 4. On the other hand, we can see from Tab.[2]that Acj(L) = u>(L) — 2 
tends to increase in magnitude at L > 13. We have illustrated this systematic and 
smooth deviation in Fig. El The only reasonable explanation of this behavior is that 
the expansion (|21j) necessarily contains the exponent 2 and, likely, also the exponent 
4, and at the same time it contains also a correction of a very small amplitude with 
to < 2. The latter explains the increase of Aui(L). Namely, the correction to scaling 
for L l ^G(L) behaves like const ■ L~ 2 [1 + (L~ 2 ) + e L 2 '^] with e < 1, which 
implies a slow crossover of the effective exponent u> e ff(L) from the values about 2 to 
the asymptotic value u. Besides, in the region where eL 2 ~ w <^ 1 holds, the effective 
exponent behaves like 

u eff (L) ~ 2 + 6iL 2 ^ + 6 2 L" 2 , (26) 

where b\ <C 1 and 62 are constants. By using the extrapolation of w e // with u' = 2 
in (|24[) and (|25[). we have compensated the effect of the correction term 62 -£ -2 - Be- 
sides, by matching the amplitude 6 in (|25|) we have compensated also the next trivial 
correction term ~ L -3 in the expansion of Lo e ff(L). It means that the extrapolated 
exponent lo{L) does not contain these expansion terms, i. e., we have 

u(L) = 2 + biL 2 "" + 5u){L) , (27) 

where 8u)(L) represents a remainder term. It includes the trivial corrections like 
L -4 , L -5 , etc., and also subleading nontrivial corrections, as well as corrections of 
order (ei 2 " w ) 2 , (eL 2 "") 3 , etc., neglected in <^). According to the latter, Eq. ijTTj) 
is meaningless in the thermodynamic limit L — > 00, but it can be used to evaluate 
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o.oo 0.08 (io/L) 5 ' 75 °- 16 

Figure 4: The exponent 2 — ui estimated from H30[) at different system sizes. From top 
to bottom (if looking on the left hand side): a = 0, 1,3.5,5.75,7.243,9.25, 12. The results 
at the optimal a value 7.243 are shown by empty circles. The dashed line indicates our 
theoretical asymptotic value 2 — lu = 1.75, whereas the dot-dashed line - that proposed 
in H7J. 



the correction-to-scaling exponent uj from the transient behavior at large, but not 
too large values of L where b\l?~^ <C 1 holds. In our example the latter condition 
is well satisfied, indeed. 

Based on (|27|). we have estimated the nontrivial correction-to-scaling exponent 
to by using the data of u>(L) in Tab[2j We have used two different ansatzs 

2-wi(JD) =ln[A£5(L)/A2(L- l)]/m[L/(L- 1)] (28) 

and 

2 - = L [Aw(L) - A£(L - 1)] /A£(L) , (29) 

as well as the linear combination of them 

ui{L) = (1-a) u) X {L) +au 2 (L) (30) 

containing a free parameter a. We have w(-L) = lo\{L) at a = and lo{L) = uo2{L) 
at a = 1. In general, the effective exponent lo{L) converges to the same result u at 
arbitrary value of a, but at some values the convergence is better. The results for 
2 — uj{L) vs L w_6 at different a values are represented in Fig. H]by a set of curves. 
In this scale the convergence to the asymptotic value would be linear (within the 
actual region where L 3> 1 and b\ L 2 ~ w <C 1 hold) for a = at the condition 
5uj(L) oc L~ 4 . We have choosen the scale of L~ 5 ' 75 , as it is consistent with our 
theoretical prediction in Sec. EJ that uj = 1/4. Nothing is changed essentially if we 
use slightly different scale as, e. g., L~ 14 / 3 consistent with the correction-to-scaling 
exponent u = 4/3 proposed in [27j. As we see from Fig. 01 all curves tend to merge 
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at our asymptotic value 2 — u = 1.75 shown by a dashed line. The optimal value of 
a is defined by the condition that the last two estimates w(17) and u;(18) agree with 
each other. It occurs at a = 7.243, and the last two points lie just on our theoretical 
line. A discussion and comparison of our results with those in published literature 
(e. g., 28 j) can be found in 



4.4 Comparison to the known exact results and estimation of nu- 
merical errors 

We have carefully checked our algorithms comparing the results with those obtained 
via a straightforward counting of all spin configurations for small lattices, as well as 
comparing the obtained values of the partition function to those calculated from the 
known exact analytical expressions. Namely, an exact expression for the partition 
function of a finite-size 2D lattice on a torus with arbitrary coupling constants 
between each pair of neighbouring spins has been reported in [301 obtained by the 
loop counting method and represented by determinants of certain transfer matrices. 
In the standard 2D Ising model with only one common coupling constant (5 these 
matrices can be diagonalized easily, using the standard techniques Besides, the 
loop counting method can be trivially extended to the cases with antiperiodic or 
mixed boundary conditions discussed in Sec. 13.31 It is necessary only to mention 
that each loop gets an additional factor —1 when it winds round the torus with 
antiperiodic boundary conditions. We consider the partition functions Z pp = Z, 
Zaa, Z ap , Z pa . In this notation the first index refers to the horizontal or x axis, 
and the second one - to the vertical or y axis of a lattice illustrated in Fig. [fli; p 
means periodic and a - antiperiodic boundary conditions. As explained above, the 
standard methods leads to the following exact expressions: 

z pp = (Qi + Q 2 + Q 3 -Qo)/2 

Zap = (Q + Ql+Q3-Q2)/2 

z pa = (Qo + Qi + Q 2 -Q 3 )/2 (31) 

Z aa = (Qo + Q2 + Q3-Ql)/2 

where Q$ is the partition function represented by the sum of the closed loops on the 
lattice, as consistent with the loop counting method in [HJ, whereas Qi, Q2, and Q3 
are modified sums with additional factors exp^Ax-in/N+Ay-iTr/L), exp(Ax-i7r/N), 
and exp(Ay -m /L), respectively, related to each change of coordinate x by Ax = ±1, 
or coordinate y by Ay = ±1 when making a loop. The standard manipulations |31| 
yield 

Q i = 2 NL J] [cosh 2 (2/3) - sinh(2/3) (32) 

Qx,q y 

■lMl/2 



X COS 



+ cos 



% + {k,\ + $i,3) y 



where the wave vectors q x = (2n/N) ■ n and q y = (2n/L) ■ I run over all the values 
corresponding to n = 0, 1, 2, . . . , N — 1 and £ = 0, 1, 2, . . . , L — 1. Eq. (|32|) represents 
an analytic extension from small (5 region j^Uj. The correct sign of square roots 
is defined by this condition, and all Qi are positive except for Qq, which vanishes 
at (3 = f3 c and becomes negative at (3 > /3 C . In the case of the periodic boundary 
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conditions, each loop of Qo has the sign (^—\^ m + ab + a + b |3Qj . where m is the number 
of intersections, a is the number of windings around the torus in x direction, and b 
- in y direction. The correct result for Z pp is obtained if each of the loops has the 
sign (— l) m . In all other cases, similar relations are found easily, taking into account 
the above defined additional factors. Eqs. (|31|) are then obtained by finding such a 
linear combination of quantities Qi which ensures the correct weight for each kind 
of loops. 

All our tests provided a perfect agreement between the obtained values of the 
Greens functions G(r) (a comparison between straightforward calculations and our 
algorithms), as well as between partition functions for different boundary conditions 
(a comparison between our algorithms and Eq. I)31jl ). The relative discrepancies were 
extremely small (e. g., 10 -15 ), obviously, due to the purely numerical inaccuracy. 

We have used the double-precision FORTRAN programs. The numerical errors 
in Tab. 121 have been estimated by repeating some calculations with twice larger 
number of digits (REAL*16 option). Thus, the errors in the G(L) values for L = 10 
to L = 17 are 4.7 • 10~ 17 , 4.06 • 10" 16 , -3.52 • 10~ 16 , -5.65 • 10~ 16 , 1.03 • 10~ 15 , 
1.41 • 10 -15 , —1.71 • 10 -16 , and 3.09 • 10~ 16 . To eliminate the summation error for 
the largest lattice L = 20, we have split the summation over the configurations of 
the first row in several parts in such a way that a relatively small part, including 
only the first 10 000 configurations from the total number of 52 487 nonequivalent 
ones, gives the main contribution to Z and Z'. The same trick with splitting in two 
approximately equal parts has been used at L = 19. As a result, the numerical errors 
at L = 18, 19, 20 are not much larger than the above listed values for 10 < L < 17. 
Hence, the resulting numerical errors in Fig.0]do not much exceed 0.03 in the middle 
part around 2 — uj ~ 1.75. In Fig. the errors are practically not seen. 

5 Generation of pseudo— random numbers 

We have found that some of simulated quantities like specific heat of 3D Ising model 
near criticality are rather sensitive to the quality of pseudo-random numbers. The 
linear congruatial generators providing the sequence 

In+i = (aln + c)mod m (33) 

of integer numbers I n is a convenient choice. We have used in previous section 
the generator of 33 with a = 843314861, c = 453816693, and m = 2 31 . The 
G05CAF generator of NAG library with a = 13 13 , c = 0, and m = 2 59 (generating 
odd integers) has been extensively used in [33] . We have compared the results of 
both generators for 3D Ising model, simulated by the Wolff's cluster algorithm 36 , 
and have found a disagreement by almost 1.8% in the maximal value of Cy at the 
system size L = 48. Application of the standard shuffling scheme ( |38| p. 391) with 
the length of the shuffling box (string) N = 128 appears to be not helpful to remove 
the discrepancy. The problem is that the standard shuffling scheme, where the 
numbers created by the original generator are put in the shuffling box and picked 
up from it with random delays in about N steps, effectively removes the short- 
range correlations between the pseudo-random numbers, but nevertheless it does 
not essentially change the block-averages {I n )k = Y^j=n~ l over ^ subsequent 
steps if k » N. It means that such a shuffling is unable to modify the low-frequency 
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tail of the Fourier spectrum of the sequence I n to make it more consistent with white 
noise (an ideal case). The numbers I n repeat cyclically and the block-averages over 
the cycle do not fluctuate at all in contradiction with truly random behavior. To 
solve the problem, we have made a second shuffling as follows. We have split the 
whole cycle of length m of the actual generator with m = 2 31 in 2 20 se gments 
each consisting of 2048 numbers. Starting with 0, we have recorded the beginning 
numbers of each segment. It allows to restart the generator from the beginning of 
any segment. The last pseudo-random number generated by our shuffling scheme is 
used to choose the next number from the shuffling box, exactly as in the standard 
scheme. In addition, we have used the last but one number to choose at random 
a new segment after each 2048 steps. This double-shuffling scheme mimics the 
true fluctuations of the block-averages even at k 3> m and has an extremely long 
(comparable with (^o) (2 20 )! steps at N = 2 20 ) cycle. We have used a very large 

shuffling box with N = 2 20 to make the shuffling more effective. As a consequence, 
we have reached a perfect agreement with the results of G05CAF generator, which 
has a rather long cycle even without shuffling. 

A hidden problem is the existence of certain long-range correlations in the se- 
quence I n of the original generator of [HE]- Namely, pseudo-random numbers of a 
subset, composed by picking up each 2 fc -th element of the original sequence, appear 
to be rather strongly correlated for k > 20. It is observed explicitely by plotting 
such a subsequence I* vs n, particularly, if the first element is choosen I* = 0. 
These correlations reduce the effectiveness of our second shuffling. Correlations of 
this kind, although well masked, exist also in the sequence of G05CAF generator. 
Namely, if we choose I* = 1 and k = 25 and generate the coordinates (x, y) by 
means of this subset, then we observe that the x — y plane is filled by the generated 
points in a non-random way. The origin of these correlations, obviously, is the choice 
of modulo parameter m as a power of 2. It, evidently, is the reason for systemat- 
ical errors in some applications with Swendsen-Wang algorithm discussed in |37| . 
A promising alternative, therefore, is to use the well known Lewis generator 38 , 
where m = 2 31 — 1 is a prime number, a = 7 5 , and c = 0, as the original generator 
of our double-shuffling scheme. (This generator has been tested in jHS] and, even 
without any shuffling, it gave good results for the energy and specific heat of 2D 
Ising model on 16 x 16 lattice simulated by Wolff's cluster algorithm.) As before, 
the cycle is split in 2 20 segments. However, the first segment now starts with 1. 
Besides, the first and the last segments contain only 2047 elements instead of 2048. 
After all numbers of the previous segment are exhausted, a new segment is choosen 
as follows: if the last but one random number of our shuffling scheme is I, then we 
choose the fc-th segment, where k = 1 + [//2048]. Since we never have / = or 
/ = m, it ensures that each segment is choosen with the probability proportional to 
its length. We have used the shuffling box of length N = 10 6 for this scheme. 

From the theoretical point of view, the latter scheme could provide the best 
pseudo-random numbers. The test simulations we made in 2D Ising model showed 
that G05CAF generator, as well as both shuffling schemes provide very accurate 
results, which indicates that the actually discussed long-range correlations do not 
have a remarkable effect in our application. We have simulated by the Wolff's 
algorithm the mean energy (e), specific heat Cy, as well as its derivatives C' v = 
dCv/d/3 and C v = d 2 Cyjd0 1 for 2D Ising model at the critical point and have 
compared the results with those extracted from exact formulae ()31j) and (|32[). The 
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Figure 5: A test estimation of the critical exponent (3 in 2D Ising model by the method 
of effective exponents (left) and by measuring the slope of the InM vs hit plot (right). 
The tiny-dashed and the solid lines in the left-hand-side picture show the linear and the 
quadratic extrapolations yielding j3 ~ 0.1244 and (3 ~ 0.12496, respectively. The horizontal 
dashed line shows the exact value 1/8. The slope of the linear 6-point fit within t £ 
[t m in',t max ) in the right-hand-side picture gives f3 ~ 0.1144 ~ /3 e //(t*), where the values 
0.1144 and t* = \/t m i n t max are indicated in the other picture by dot-dot-dashed lines. 

test simulations consisting of 4.8- 10 8 and 2.4- 10 7 cluster-algorithm steps have been 
made for the lattice sizes L = 48 and L = 256, respectively. The whole simulation 
has been split in 24 blocks to estimate the statistical averages and standard errors 
(a) from the last 20 blocks. The simulation with the generator of jSSl has revealed 
systematical errors of about 10a for the specific heat and its derivatives at L = 48. 
The values provided by the G05CAF generator and our two shuffling schemes agreed 
with the exact ones within the errors about one a. The most serious deviation of 
2.37a has been observed for Cy in the case of L = 48 simulated by our first shuffling 
scheme. At L = 48, one standard deviation a corresponded to ~ 0.0009% relative 
error for (e), ~ 0.02% error for Cy, ~ 0.2% error for Cy, and ~ 0.35% error 
for C v . At L = 256 these errors were ~ 0.0012%, ~ 0.12%, ~ 3%, and ~ 4%, 
respectively. Furthermore we have used the latter three generators in simulations of 
3D Ising model and verification of the simulated values by performing some of the 
simulations twice with different generators. 

6 Test estimations of the critical exponent f3 in 2D Ising 
model 

Based on the well known exact mag netization data M = (1 - [sinh 2K]- 4 f of 2D 
Ising model, we have tested the known method of effective exponents |40M41| exten- 
sively used in our paper. Here (3 = 1/8 is the magnetization exponent and K is the 
coupling constant, denoted in this way exceptionally in Sees. El to El 

The effective critical exponent (3 e ff(t) = ln[M (at) /M(t/ a)]/ (2 In a) with a = 
2 1//4 , calculated from the magnetization data M(t) within the range of the reduced 
temperatures t = (K/K c ) — 1 € [0.02; 0.08] (where K c is the critical coupling) is 
shown in Fig. El (left) by solid circles. Omitting two largest t values, the linear 
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least-squares approximation of /3 e //(t) (tiny dashed line) gives j3 ~ 0.1244, and the 
quadratic fit of all points (solid line) yields (3 ~ 0.12496 in close agreement with 
the exact value 0.125. For comparison, the most popular method of estimation of 
critical exponents by simply measuring the slope of a log-log plot, as illustrated in 
the right-hand-side picture (Fig. |5J), yields a relatively poor result (3 ~ 0.1144 in 
spite of the fact that the actually used piece of the log-log plot (6 smallest t values) 
looks very linear. Up to now we have discussed the exact data only. In the case of 
Monte Carlo data with the statistical errors, say, about the symbol size, we would be 
unable to detect the very small deviations from linearity and could easily get a very 
good, but nevertheless misleading linear fit. In other words, such a simple measuring 
of critical exponent is unreliable since there is clearly a danger to get an uncontrolled 
systematical error. Such a measurement within t £ [t m i n ', t ma x\ practically yields the 
mean slope of the log-log plot within this interval, which is nothing but an effective 
exponent. It corresponds just to one point on the (3 e ff{t) plot, i. e., to (5 e ff{t*) at 
t* = \/tmi n tmax, as indicated in Fig. [5] (left) by dot-dot-dashed lines. The method 
of effective exponents has been designed to control the systematical errors of such 
simple measurements and eliminate them by a suitable extrapolation. Inclusion of 
corrections to scaling directly in the ansatz for row data (in this case M) is another 
way to eliminate the systematical errors. However, we greatly prefer the method 
of effective exponents since it can be well controlled visually. This method is also 
sensitive enough to distinguish between a power-like and a logarithmic singularity 
of specific heat, as discussed in Sec. El 

There are no essential problems reported in literature, as regards the MC estima- 
tion of the critical exponent (3 in 2D Ising model. It is because the simulations can be 
done easily much closer to the critical point than in our test example. However, the 
problem remains in 3D case. The systematical errors in the measured (3 values are 
caused by the corrections to scaling, so that t rather than t is an essential param- 
eter. For the smallest reduced temperatures t > 0.0005 considered in the published 
literature |l2] we have t e > 0.022 with the RG value of the correction-to-scaling 
exponent 6 ~ 0.5, and t 9 > 0.079 with our (GFD) value 9 = 1/3. It means that the 
systematical errors of the simple (naive) measurements of (3 can be even larger than 
in our 2D test example with t s = t > 0.02. 

7 Estimation of the critical coupling of 3D Ising model 

Here we discuss the critical coupling K c of 3D Ising model, which is relevant to 
our estimations of the critical exponent (3 in Sec. |H1 The most accurate MC values 
reported in literature are K c = 0.22165459(10) gSj, K c ~ 0.2216545 jHj, and K c = 
0.2216544(3) One of the recent estimates of HT series expansion is K c = 

0.221654(1) We have estimated the critical coupling from our MC results 

for the pseudocritical couplings K C (L) which correspond to U = 1.6, where U = 
(M 4 )/(M 2 ) 2 . Note that 1 - U/3 is the Binder cumulant which may have the values 
from (at high temperatures) to 2/3 (at low temperatures). In the thermodynamic 
limit L — > oo it changes jump-likely at K = K c , so that lim K C (L) = K c holds for 

L— >oo 

any given U within 1 < U < 3. 

The values K c (48) = 0.22164540(118), K c (64) = 0.22165095(153), K c (96) = 
0.221653069(734), K c (128) = 0.221653945(453), K c (192) = 0.221654550(316), 
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K c (256) = 0.221654755(163), and K c (384) = 0.221654532(109) have been obtained 
by an iterative method similar to that described in Sec. |HJ 

The data suggest that the pseudocritical coupling has a maximum at L ~ 256. 
Since dU/dK is negative, such a qualitative behavior is expected in view of the 
known results [12], according to which the universal value of U at K = K c and 
L — ► oo is slightly larger than 1.6 and, therefore, K C (L) should approach K c from 
above. According to the finite-size scaling theory, K C (L)-K C oc L~ x l v [1 + (L~ w )] 
holds at large L, where to = 9 /v. We have found that the data within a wide range 
of sizes can be well described by the Pade approximation rather than by a simple 
ansatz with the correction-to-scaling term. Namely, the formula 

K C (L)*K C + C-V \ + + a £," , (34) 

where C = L/Lq and Lq corresponds to the maximum of K C (L) plot, well fits the 
data within the whole range of sizes L € [48; 384]. The location of the maximum is 
the only characteristic length measure for the K C {L) plot, which should transform 
to the correct asymptotic form somewhat beyond this maximum. It motivates our 
specific choice of Lq, which otherwise is not well defined as a fitting parameter. 
Fortunately, the results remain practically unchanged if we take, say, twice smaller 
or twice larger value of Lq. 

Assuming our (GFD) exponents v = 2/3 and to = 0.5 (correction-to-scaling 
exponent for the magnetization), a fitting of all data points to Q34JI yields K c = 
0.22165407(29) with the goodness of fit [H] Q = 0.897. To eliminate the systematical 
errors, we have discarded the two smallest sizes, which yields K c = 0.22165386(51) 
with Q = 0.797. It agrees within the error bars with the value 0.2216544(3) of [22] 
and the value 0.221654(1) of |45| . Our value is provided by the fit within L £ [96; 384] 
at Lq = 234. It is shifted up (down) only by 5.5 • 10 -8 (6.1 • 10 -8 ) at a twice smaller 
(larger) Lq. It is also not very sensitive to the choice of the exponents v and u. 
Assuming the RG exponents v = 0.63 and uj = 0.8, we obtain K c = 0.22165395(46) 
with Q = 0.795. Further we have used both our values of K c and those reported in 
literature in various estimations of the critical exponent (5. 

8 Estimation of the critical exponent (3 from the mag- 
netization data in 3D Ising model 

Based on the well known scaling relation 2(3 = dv — 7, we find from ([2]) and (jHJ), 
where j = and m = 3 holds at n = 1, the GFD value (5 = 3/8 of the magnetization 
exponent /3 for 3D Ising model. This value is remarkably larger than the usually 
believed ones about 0.326 [T2|. We suppose that the asymptotic exponent (3 not 
only for the Ising model, but also for the Heisenberg model is larger than provided 
by approximate RG theories and known numerical estimates. In poly crystalline 
Ni (n = 3), the increase of the effective exponent j3 e ff well above the RG value 
0.3662 [17] has been established experimentally in |48] . where the authors have found 
the asymptotic estimate [3 = 0.390(3). This value clearly disagree within error bars 
with the RG prediction, but agree with our value 11/28 = 0.3928. . . predicted for 
the n = 3 case (m = 3, j = 2). Also the critical exponent 7 measured in most of 
experiments on Ni and Fe ranges from 1.28 to 1.35 (see |48[ I49 | IBTi] and references 



20 



therein), and only some experimentators have obtained a larger value about 1.41 - 
the only value cited in |12j . One of the best experimental methods is to use the 
Kouvel-Fisher plot, since T c and 7 are determined simultaneously with no fitting 
parameters 0§1- This method yields the value 7 = 1.35 |491 151j which is believed 
among (some) experimentators to be the asymptotic exponent (see the references 
in |49j). Our prediction 7 = 19/14 ~ 1.357 is remarkably consistent with this value. 

Let us now return to the Ising model. The spontaneous magnetization M of 3D 
Ising model has been considered in |521 1531 IS"4l I42j . An empirical formula 

M(t) = t - 32694109 (l.6919045 - 0.34357731 f - 50842026 - 0.42 5 7 2 3 66 tj (35) 

with t = 1 — 0.2216544/7^ has been found in |42| which fits the simulated at 
three linear sizes L = 64, 128, 256 of the lattice and extrapolated to the ther- 
modynamic limit data of (| M |) within the range of the reduced temperatures 
t ~ t = (K/K c ) — 1 > 0.0005. We have made an approximate estimatiom of the 
spontaneos magnetization M(t) in the thermodynamic limit from our simulated val- 
ues of (| M |) at L = 200 to compare the results with (|35|) . Besides, we have made 
accurate MC simulations by the Wolff's algorithm pHU at t > 0.000086 for system 
sizes up to L = 410 to verify the critical exponent (3 ~ 0.3269 proposed by (|35|) . 
We have performed the simulations at certain values of coupling constants Ki, 

n 224 — K* 

Ki = K + » "2<^<14 (36) 

(rounded to 9 digits after the decimal point), where K* = 0.2216545 is the critical 
coupling estimated in We have choosen Kq = 0.224 to compare our results 

with those reported in jH3H^. Each next Ki value is v2 times closer to K* than 
the previous one. 

Our MC simulations have been split typically in 51 blocks (bins) to calculate 
the average value an standard deviation a from the last 50 bins. As a test, we 
have checked that a splitting in twice larger blocks provides consistent values of a, 
i. e., the blocks were large enough to justify our treatment of the block-averages 
as statistically independent quantities. In most of the cases one bin included J = 
120000 cluster updates, which corresponds to J(M 2 ) complete updates of the whole 
system or sweeps, as consistent with the known improved cluster estimator (M 2 ) = 
L~ d (c) (HE], where (c) is the average cluster size. Somewhat shorter simulations have 
been performed for our estimations of (| M |) in Tab.|3 The results at our smallest 
K value = 0.221672824 have been obtained by an averaging over four runs, i. e., 
4 x 50 x 60000 cluster updates. Note that the auto-correlation time of the Wolff's 
algorithm at criticality is only few (2 or 3) sweeps [HE]- The MC measurements 
were made frequently with respect to this auto-correlation time, i. e., the fraction of 
moved spins between the measurements was about 0.15 or smaller. In all cases we 
have discarded no less than 300 sweeps from the beginning of the simulation. We 
have verified that the system has been equilibrated well enough by comparing the 
estimates from separate smaller parts of the whole simulation. 

Our (| M |) data for K^2 < K < K3 are listed in Tab. El The second shuffling 
scheme described in Sec. El has been used as a source of pseudo-random numbers 
for this simulation. The values of (| M |) are only slightly varied with L, and those 
at the largest two sizes L = 120 and L = 200 agree or almost agree within the error 
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Table 3: The simulated values of (| M |) at four different system sizes L in the range of 
coupling constants K 3 < K < K- 2 - 





L = 60 


L = 80 


L = 120 


L = 200 


K = K_ 2 
K = K_i 
K = K 
K = Ki 
K = K 2 
K = K 3 


0.460754(205) 
0.414874(203) 
0.372813(339) 
0.334641(378) 
0.299917(295) 
0.268433(402) 


0.460695(149) 
0.414845(236) 
0.372785(240) 
0.334207(221) 
0.299268(298) 
0.268930(294) 


0.460787(155) 
0.414329(162) 
0.372939(212) 
0.334246(211) 
0.299473(211) 
0.268807(286) 


0.460493(96) 

0.414503(141) 

0.372389(142) 

0.334327(122) 

0.299621(163) 

0.268723(200) 



bars. According to |54j . the latter size more than 22 times exceeds the correlation 
length £ at K > K$, which is quite enough to estimate the thermodynamic limit. 
Our data at L = 200 are in a reasonable agreement with the corresponding values 
0.460435, 0.414490, 0.372471, 0.334258, 0.299652, and 0.268412 given by (JHSJ) at 
K = K_ 2 to K = K z . Our result at K = K = 0.224, i. e. M = 0.372389(142) or 
M 2 = 0.138674(106), agree within the error bars with the M 2 value 0.138708(39) 
reported in |53| . as well as with M 2 value 0.1387488(75) obtained in |54j . 

Our estimation of the critical exponent (3 is based on the analysis of the effective 
exponent 

R m ln(M 2 (ti, L(ti))) - ln(M 2 (t 2 , L(t 2 ))) 

0eff{t) = 2(111*!- In t 2 ) ' (37) 

where hit = (ln*i +lnt 2 )/2, L(t)t u = const, and (M 2 (t,L)} is the statistically 
averaged squared magnetization at the given t and L. The effective exponent is 
the average slope of the 0.51n(M 2 ) vs hit plot within t £ [^i;i 2 ], calculated at a 
fixed scaling argument Lt u which corresponds to a certain asymptotic value of L/£ 
at t — ► 0, where £ ~ t~ v is the correlation length. For any given ratio t 2 /ti, the 
values of e ff(t) lie on a smooth analytical curve converging to the true asymptotic 
exponent at t — ► 0. Besides, the estimates obtained at slightly different values of 
i 2 /*i well coincide with each other. We have choosen L(t) = 256(t/t) u , where t is a 
reference value of the reduced temperature at K = K\2 = 0.221691148. In this case 
L(t) approximately corresponds to the minimum of (| M |) vs L plot 55.,, and the 
deviations from the thermodynamic limit are small. Two slightly different values 
for the correlation length exponent v have been used, i. e., u = 2/3 (our value) and 
v = 0.63 (RG value). 

We have made the simulations at L values close to L(t), which allow us to 
estimate (M 2 (t, L(t))) both at v = 2/3 and v = 0.63 by a linear interpolation 
of (M 2 ) vs L~ d plot. This plot is linear at L — > 00, as discussed in 52 j. At 
K = Kiz,Kn, all (three) points have been fit together to get a more reliable result. 
The simulated values are listed in Tab.0J In most of the cases the G05CAF pseudo- 
random number generator has been used, whereas some of the values, which are 
marked by an asterisk, have been simulated by our second shuffling scheme (Sec. Oy). 
The first value at K = K\2 and L = 256 has been obtained from 50 x 90000 cluster 
updates. The second one represents the result of extracted from the simulation 
of approximately the same length with the same G05CAF generator. As we see, 
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Table 4: Squared magnetization (M 2 (K,L)) at different coupling constants K and system 
sizes L. The simulated values have been obtained by using G05CAF pseudo-random number 
generator, except those marked by an asterisk, for which our second shuffling scheme (Sec. 01 
has been applied. The value of ^Sj is marked by a dagger. 



L 


(M 2 (K U L)) 


L 


(M 2 (K , L)) 


L 


(M 2 (K 1 ,L)) 

\ \ 1' // 


1 9 




-LU 


U. 14:00001 l^U 1 


9D 


U. J. ±c/Q£0\ 104: 1 


13 


0.182897(143) 


18 


0.144657(144) 


21 


0.117712(98) 


14 


0.180291(156) 


19 


0.143549(125) 


23 


0.116133(116) 


15 


0.177886(154) 






24 


0.115459(127) 


L 


(M 2 (K 2 ,L)) 


L 


(M 2 (K 3 ,L)) 


L 


(M 2 (K 4 ,L)) 


95 


fl 095551 ('791 


39 


D 076381 ('99") 


4D 


D 060939(1 05"l 


26 


0.095077(113) 


32 


0.076290(109)* 


41 


0.060583(87) 


28 


0.093522(99) 


35 


0.074865(99) 


44 


0.059919(68) 


29 


0.092933(107) 


36 


0.074543(103) 


45 


0.059726(92) 


L 


(M 2 (^ 5 ,L)) 


L 


<M 2 (K 6 ,L)> 


L 


(M 2 (K 7 ,L)) 


50 


0.048762(79) 


64 


0.038874(80) 


80 


0.031042(79) 


51 


0.048476(89) 


64 


0.038872(69)* 


83 


0.030865(68) 


55 


0.047911(79) 


69 


0.038338(83) 


86 


0.030618(69) 


56 


0.047747(75) 


70 


0.038252(79) 






L 


<M 2 (^ 8 ,L)> 


L 


(M 2 (K 9 ,L)) 


L 


(M 2 (K 10 ,L)) 


101 


0.024673(66) 


128 


0.019726(73) 


161 


0.015670(49) 


104 


0.024538(60) 


128 


0.019735(50)* 


166 


0.015575(46) 


107 


0.024466(63) 


133 


0.019596(59) 






L 


{M 2 {K^L)) 


L 


(M 2 (K 12 ,L)) 


L 


(M 2 (K 13 ,L)) 


203 


0.012498(45) 


256 


0.009936(47) 


315 


0.007950(32) 


206 


0.012464(62) 


256 


0.009862(49) t 


320 


0.007911(43) 






256 


0.009976(87)* 


320 


0.007943(39)* 










325 


0.007901(38) 


L 


(M 2 (K 14 ,L)) 


L 


(M 2 (K U ,L)) 


L 


(M 2 (K 14 ,L)> 


390 


0.006344(25) 


400 


0.006317(27) 


410 


0.006275(25) 
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Figure 6: The effective critical exponent /3 e // vs t e with 9 = 0.5 (left) and = 1/3 
(right). The values calculated from with v = 0.63, K c = 0.22165395 (left) and v = 2/3, 
K c = 0.22165386 (right) are shown by solid circles (averaged over K = Ki+j, Ki+^-j with 
j = 0, 1, 2) and empty circles (averaged over K = Ki+j, Ki+4-j with j = 0, 1). Solid lines 
represent the cubic fits, whereas the dashed line (left) shows the linear 12-point fit. The 
dot-dot-dashed line shows the RG prediction j3 — 0.3258 fT7) . 

different simulations confirm each other within the error bars. We have made a 
weighted averaging (with the weights proportional to the simulation length which 
is roughly oc l/c 2 ) of the over taping simulation results to obtain statistically more 
reliable values for our further analysis. 

We have calculated from 1|37|) and plotted in Fig. |f)]the effective critical exponent 
(3 e ft as a function of t with the RG exponents v = 0.63 and 6 = 0.5 (left), as well 
as with the exponents of GFD theory v = 2/3 and 9 = 1/3 (right). The corre- 
sponding self consistent estimates of the critical coupling K c = 0.22165395(46) and 
K c = 0.22165386(51) have been used, obtained in Sec. [7] from the Binder cumulant 
data within L £ [96; 384]. The results of a weighted averaging over the estimates 
obtained at K = Ki+j, Ki+^-j with j = 0, 1, 2 are shown by solid circles, whereas the 
values averaged over K = Ki+j, Ki+^-j with j = 0, 1 are depicted by empty circles. 
The f3 e f f values at K = Kj , Kj+i have been taken with the weights oc £ 2 , which 
approximately minimize the resulting statistical errors, taking into accont that the 
individual errors are roughly proportional to 

In both cases (left and right) the plot of the effective exponent has an inflection 
point and is well described by a cubic curve, although the last 12 points (smallest 
t values) can be well fit by a straight line, too. The linear dashed-line fit with the 
RG exponents (left) yields (3 = 0.3347(24) at a fixed K c = 0.22165395. Taking 
into account the uncertainty in K c , we have (5 = 0.3347(52). This value reduces 
to 0.3302(37) if we take K c = 0.2216544(3) estimated in @2j. Nevertheless, in 
both cases it is slightly larger than the RG value 0.3258 07j (dot-dashed line), 
supported by the high temperature (HT) series expansion 45 , and also a bit larger 
than the value of [121 P = 0.3269(6) [ansatz (|35|) ]. The linear fit, in fact, takes into 
account the leading correction to scaling only. The cubic fit at K c = 0.22165395(46) 
(solid curve), which includes also two sub-leading corrections, tends to deviate up 
to P = 0.3385(73) (0.3385(37) at a fixed K c ) in a worse agreement with the RG 
prediction. Assuming the critical coupling K c = 0.2216544(3) of |32j, the cubic fit 
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gives (3 = 0.3324(52). It is worthy to mention that the fits supporting the RG value 
with a striking accuracy can be produced easily without simulations so close to the 
critical point. For instance, omitting the 5 smallest K values (10 points on the 
f3 e ff plot), which roughly corresponds to the simulation range t > 0.0005, the linear 
10-point fit yields (3 = 0.3262(16) at K c = 0.2216544(3) (the K c value of 42 j) and 
(3 = 0.3259(15) at K c = 0.22165459(10) (the K c value of @3]). 

Taking into account that the j3 e ff data are correlated, the statistical errors of 

1 /9 

the extrapolated values have been estimated as 8f) > where Si is the partial 
error due to the uncertainty in the i-th value of (M 2 ). 

A self consistent extrapolation within the GFD theory is illustrated in the right 
hand side picture (Fig. EJ. The cubic fit gives [3 = 0.366(16) in agreement with the 
expected exact result 0.375. Unfortunately, there is still a very large extrapolation 
gap, so that we cannot make too serious conclusions herefrom. It is necessary to go 
even much closer to the critical point in this case with simultaneous reduction of the 
error in the estimated K c value. 

9 Estimation of the singularity of specific heat in 3D 
Ising model from the finite— size scaling of MC data 

It is commonly believed [Uj that the specific heat Cy of 3D Ising model on an in- 
finitely large lattice has a power-like singularity, i. e., Cy ~ t~ a at t — * 0. According 
to the finite-size scaling theory, it would mean that 

Cy ~ t~ a f (L l / y t) = L a ' u f (L^tj (38) 

holds at small t in the finite-size scaling region t ~ L~ x l v . Here v is the exponent of 
correlation length, whereas f(z) and f{z) = z~ a f{z) are the scaling functions. The 
maximum of the Cy vs t plot is located at a certain value of the scaling argument 
z = I?-I v t at L — > oo, which would mean that the maximum values scale as 

C^ X (L) oc L a ' v at L -> oo . (39) 

An estimation of the exponent a/u from the slope of the log-log plot then gives us 
the effective exponent (a/v) e ff = <91nC™ ax (L)/dlnL, which is varied due to the 
corrections to scaling like 

(a/v) eff ~ (a/u) + const ■ , (40) 

where u = 9/v is the correction-to-scaling exponent. Note that specific heat con- 
tains also an analytic background contibution which influences this behaviour. We 
have found, however, that this influence is very small in the actually considered 
range of sizes if the constant background term is of order one, as expected from 
physically-intuitive considerations. 

We allow a possibility that the specific heat has a logarithmic singularity, as 
consistent with [201 an d (Ej- It means that Eq. (j3*H|) is replaced with 

C v ~\nt- f (41) 
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Figure 7: The effective exponent (a/v) e ff of 2D Ising model depending on x(L) = 
l/ln(L/Lo) with Lq = 0.572 (solid circles) and x(L) = (empty circles). The solid 

straight line represents the approximation (pt/v) e ff = l/ln(L/Lo). 



and (HOJ) - with 

{a i v) *>> k Rik) ■ (42> 

where Lo is a constant length scale. Although one believes usually that the singu- 
larity of Cy is power-like with a — 0.11, no strong numerical evidences exist which 
could rule out the logarithmic singularity (|41|). The problem is that hit behaves 
almost like a weak power of t with the effective exponent 1/lnt, e. g., like t~ 
at t ~ 10~ 4 . Moreover, below we will show that the finite-size scaling of the maxi- 
mal values of Cy is even very well consistent with (1421) in favour of the logarithmic 
singularity. 

We have tested our method in 2D Ising model, where the logarithmic singularity 
of specific heat is well known. We have considered the effective exponent 

(a/v) eff {L) = In [C^(2L)/Cy™(L/2)]/\n4 , (43) 

which is defined by finite differences of the log-log plot. It coincides with 
d In C™ ax (L)/<9 In L at large L where the log-log plot is almost linear. Based on 
exact data, extracted from (|31j) and (|3*2*|) . we have found that the effective expo- 
nent {a/v) e ff{L) within L G [48; 512] is fairly well described by ansatz (|42|) with 
Lq = 0.572 and remarkably worse described by ansatz (|40[). It is evident from Fig.0 
where {a/u) e ff vs x(L) = 1/ In (L/Lq) plot (solid circles) well coincides with the 
theoretical straight line and is much more linear than the (a/v) e ff vs x(L) = 10L -1 
plot (empty circles). Thus, our method allows to distinguish between the logarith- 
mic singularity and a powerlike singularity including correction to scaling of the kind 
oc where u = 1 holds in 2D case. 

The specific heat as well as its derivatives with respect to the coupling constant 
can be calculated easily from the Boltzmann's statistics. Thus, omitting an irrelevant 
prefactor, the specific heat is given by 

Cy = N ((e 2 ) - (e) 2 ) , (44) 
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and the derivatives of (|44|) are 




N 2 ( 3 ( e >( e 2>-( e »)-2( e > 8 ) 

A^ 3 (l2(e) 2 (e 2 ) - 3{e 2 ) 2 - 4{e){e 3 ) - 6(e) 4 + (e 4 )) 



(45) 



(46) 



where N = L 3 is the total number of spins and e is the energy per spin. The 
maximum of specific heat is located at a pseudocritical coupling (3 C which is defined 
by the condition dCy /d(5 = 0. It can be found by the Newton's iterations 



where (3 C denotes the m-th approximation of f3 c , and the derivatives are calculated 
from (jlHJi and fflty at B = $ n) . 

We have used the Wolff's algorithm (in 3D case) to estimate these derivatives in 
each iteration consisting of either 5 • 10 5 (at L < 48) or 10 6 MC steps. Besides, the 
first iteration has been used only for equilibration of the system retaining the initial 
estimate of /3 C . After few iterations (5 {Cy) reaches /3 C (C™ 1 ") within the statistical 
error and further fluctuates arround this value. In principle, the fluctuation ampli- 
tude can be reduced to an arbitrarily small value by increasing the number of MC 
steps in one iteration. 

An obvious advantage of this iterative method is that the maximal value of Cy 
can be evaluated in one simulation without any intermediate analysis. Omitting 
first 5 iterations, the mean values and the standard deviations have been evaluated 
by jakknife method [HE] from the rest 19 iterations, except the largest sizes 64 < 
L < 128, where only 4 iterations (with twice larger number of MC steps) have been 
discarded and 11 iterations have been used for the estimations. In the most of the 
cases the first shuffling scheme discussed in Sec. has been used as a source of 
pseudo-random numbers, and the simulated values have been verified by repeating 
the simulations at L = 24, 48, and 96 with the G05CAF generator. The perfect 
agreement confirms our results. 

We have averaged the values of C™ ax over both simulations at L = 24, 48, and 
96 to reduce the statistical errors in the estimated effective critical exponent (|43JI . 
Our results are summarized in Tab. El The effective exponent (a/v) e ff(L) within 
L 6 [12; 64] is rather well approximated by (|42[) with Lq = 1.258, as shown in Fig. |H] 
by solid circles and linear least-squares fit in the scale of x(L) = l/ln(L/Lo). This 
is a remarkable agreement, taking into account that ansatz 1)42(1 contains only one 
adjustable parameter Lq. We have tested also ansatz ([40)1 . where the exponents 
a/v = 0.173 and u = 0.8 have been taken from the perturbative RG theory |4JJ. 
In this case the agreement with the data is worse, i. e., the mean squared deviation 
is 5.14 times larger, as shown in Fig. [S] by empty circles and linear least-squares fit 
in the scale of x(L) = 2L~ 08 . It can be well seen also when comparing the linear 
fit of circles with the evidently better dashed-line fit. The latter represents the 
lower straight line in the scale of 2L~°' & and shows the expected behavior of empty 
circles at L — > oo if (|42jl is the correct ansatz. Note that the main deviations of 
the data points from the fitted lines in Fig. |H] are not of statistical character, since 
the statistical errors are remarkably smaller than the symbol size except only the 



£(n+l) = ~ {n) 



dCy/d/3 
d 2 C v /dp : 



;2 ' 



(47) 
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Table 5: The MC estimates of the maximal values of the specific heat C™ ax , the pseudocrit- 
ical couplings (3 C , and the effective critical exponents {a/v) e ff depending on the system size 
L. The marked values have been simulated by G05CAF pseudo-random number generator. 
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/imax 

u y 




\ u / "Jeff 


3 


16.4445(61) 


0.233595(31) 




4 


21.532(10) 


0.234207(26) 




6 


28.908(27) 


0.231090(33) 


0.66742(89) 


8 


34.155(28) 


0.228561(23) 


0.5552(14) 


12 


41.481(49) 


0.225771(20) 


0.4464(12) 


16 


46.491(85) 


0.224436(13) 


0.3894(19) 


24 


53.71(10) 


0.223207(13) 


0.3333(19) 


24 


53.65(10)* 


0.223193(16)* 




32 


58.60(15) 


0.2226726(96) 


0.3096(20) 


48 


65.82(25) 


0.2222037(97) 


0.2784(23) 


48 


65.86(18)* 


0.2221990(87)* 




64 


71.41(15) 


0.2220053(85) 


0.2593(69) 


96 


78.91(23) 


0.2218379(56) 




96 


79.01(39)* 


0.2218387(55)* 




128 


83.95(78) 


0.2217704(78) 





(a/v) eff 




0.1 0.2 0.3 0.4 „ n.5 



Figure 8: The effective exponent (a./v) e ff of 3D Ising model depending on x(L) = 
1/ ln(L/1.258) (solid circles) and x(L) = 2L~ - 8 (empty circles). The lower solid-line fit 
represents the approximation (a/v) e ff = l/\n(L/L ) with L = 1.258, whereas the upper 
one corresponds to (a/v) e ff = ipc/v) + 2.044L~" with the RG exponents ajv = 0.173 and 
uj = 0.8. The dashed line shows the logarithmic approximation in the scale of 2L~ as . Sta- 
tistical errors are within the symbol size. The mean squared deviations from the lower and 
upper solid-line fits are 1.6 • 10~ 5 and 8.22 • 10~ 5 , respectively. 
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largest L value, where the error is about the symbol size. Due to this reason we 
have used simple least-squares approximations, minimizing the sum of not weighted 
squared deviations. These deviations show just the error of the ansatz used and, 
thus, indicate that (|42|) is a better approximation than (|40j) with fixed exponents 
a/u = 0.173 and uj = 0.8 within the actual range of sizes, at least. Moreover, since 
L~ 8 has reached already a rather small value 0.036 and, therefore, the second-order 
correction ~ L~ L6 should be very small, the observed deviations can be explained 
easily assuming ()42|) rather than ()40|) with the RG exponents. According to our 
recent findings (based on unpublished simulation data for L < 192) the discrepancy 
with the RG value ajv = 0.173 can be explained by the existence of a negative back- 
ground contribution Co = —25(3) to Cy which, however, seems to be unphysically 
large in magnitude. 

It is quite possible that the true value of a = 2 — du is remarkably smaller 
than the RG value 0.11, as consistent with our recent results for the exponent \jv. 
We have estimated {l/v) e ff = y e ff(L) from the derivatives of the Binder cumulant 
U'(L) = dU(L)/d/3 (see Sec. EJ) at two system sizes L/2 and 2L. They should 
scale like U'(L) ~ L x l v at a pseudocritical coupling corresponding to U = const. 
Our data U'(W) = -175.34(0.17), U'(32) = -526.65(0.83), J7'(48) = -1007.8(2.1), 
[/'(64) = -1586.8(3.7), U'{96) = -3022.8(8.9), *7'(128) = -4773.6(16.3), J7'(192) = 
-9023.4(48.7), £/'(256) = -14284(88), and £/'(384) = -27039(171) for U = 1.6 yield 
y e //(32) = 1.5889(18), y e// (64) = 1.5901(27), y c// (96) = 1.5812(41), y e// (128) = 
1.5851(48), and y e //(192) = 1.5805(50). In analogy to the plots in Fig. El one may 
expect an accelerated further deviation from the perturbative RG value ~ 1.587. 

Thus, in spite of the conventional claims that the specific heat of 3D Ising model 
has a certain power-like critical singularity, accurately predicted by the perturbative 
RG theory, the actual very accurate MC data for C™ ax show that it is even more 
plausible that the singularity is logarithmic. 

10 Remarks about other numerical results 

There exists a large number of numerical results in the published literature not 
discussed here and in . A detailed review of these results is given in . The cited 
there papers report results which disagree with the values of the critical exponents 
we have proposed in 0. 

Particularly, the values of the perturbative RG theory are well confirmed by the 
HT series expansions [451 158j . However, we are somewhat sceptical about such a 
support of one perturbative method by another. It could well happen that the true 
reason for the agreement is the extrapolative nature of both methods, according to 
which both methods describe a transient behavior of the system far away from the 
true critical region. Really, our simulations of the magnetization data within t > 
0.0005 well confirm the RG value of the critical exponent /3, whereas the agreement 
becomes worse at t < 0.0005, where the plot of the effective exponent shows a 
remarkable inflection thus indicating that the true critical region, where the critical 
exponents can be accurately measured, is t <C 0.0005. This region, of course, cannot 
be directly accessed by the HT series expansions. Another problem is that the HT 
estimation of the critical exponent a |58j is based on a priori assumption that the 
singularity of specific heat is power-like, whereas the MC data (Sec. |HJ) suggest that 
it, very likely, is logarithmic. 
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According to the finite-size scaling theory, tL l l u is a relevant scaling argument, 
so that not too small values of the reduced temperature t are related to not too 
large system sizes L ~ t~ v . Thus, according to the idea proposed above, it is quite 
possible that the MC results for finite systems, like also the simulations at finite 
t values, appear to be in a good agreement with the conventional RG exponents 
which are valid within a certain range of t and L values well accessible in MC 
simulations. The huge number of numerical evidences in the published literature 
(see [57j) for the exponents of the perturbative RG theory certainly is a serious 
argument. Nevertheless, there is a reason to worry about the validity of these 
exponents at t — > and/or L — > oo because of the following problems. 

• We have made accurate MC simulations of the magnetization (Sec. EJ) for 
unusually large system syzes (L < 410) much closer to the critical point (t > 
0.000086 instead of t > 0.0005) than in the published literature, and have 
found that the agreement with the RG exponent (5 becomes worse in this case. 

• Our MC estimation of the exponent discussed at the end of Sec. El shows 
a good agreement with the RG value 1.587 at not too large system sizes. 
However, the agreement becomes worse when larger than L = 128 sizes are 
included. 

• A remarkable deviation of the correction-to-scaling exponent uj from the per- 
turbative RG value u ~ 0.8 has been already reported in literature [52], where 
also larger than usually system sizes L < 256 (instead of the conventional 
L < 128 [HE]) have been simulated in application to the Monte Carlo renor- 
malization group techniques, yielding uj ~ 0.7. 

• The confirmation of RG exponents by MC simulations is not unambiguous. 
There exist also examples where the simulation results are in a remarkable 
disagreement with these exponents even within the conventionally considered 
range of sizes and reduced temperatures. A particular example is the finite- 
size scaling of the maximal values of the specific heat considered in Sec. El 
We are afraid that there are also other such examples, but they are routinely 
ascribed as unbelievable and do not appear in the published literature. 

• It is indeed easy to produce evidences supporting the RG exponents, as we 
have shown in Sec. El just because these exponents describe the behavior of a 
system not too close to criticality, in the range which can be easily accessed in 
MC simulations. The problem is that the usual MC measurements yield only 
effective exponents, as shown in Sec. El which exhibit quite large variations also 
in 3D it is particularly well seen from our plots of the effective expo- 
nents. The leading correction-to-scaling term, included in the fitted ansatz, 
also does not completely solve the problem: in essence it is the same as to 
make a linear extrapolation of the effective exponent, but, e. g., the (3 e ff(t) 
plots in Fig. El are remarkably nonlinear. Therefore, only such evidences are 
really serious, which show very precisely how the effective exponents provided 
by simple estimations converge to certain asymptotic values. 

If one consider seriously a possibility that the true values of the critical exponents 
are those proposed in 0, then a question arises why the published MC estimates 
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tend to deviate greatly from these theoretical values. In our opinion, the main reason 
is that the published simulations have been made too far away from the true critical 
region (as regards both t and L), where the critical exponents can be precisely 
measured in a simple way routinely used in MC analysis. The plots of the effective 
exponents in Fig. H3 provide an evidence for this statement: as we have already 
mentioned (Sec. EJ), simple MC measurements yield just such effective exponents, 
and they are varied. One has to consider corrections to scaling, and not only the 
leading one, to get better results. However, all the existing (MC) correction-to- 
scaling analyses in the published literature rely on the RG correction-to-scaling 
exponents, therefore the disagreement with the predictions in jKj is not surprising. 

Finally, our theory provides a self consistent explanation why much smaller t 
values and/or much larger system syzes L has to be considered, as compared to the 
known simulations. It is because the correction-to-scaling exponent oj = jv = 0.5 
in our theory is remarkably smaller than that of the RG theory uj ~ 0.8, which 
implies that the decay of corrections to scaling is relatively slow. In fact, the reduced 
temperatures we have reached in our simulations of the magnetization also are still 
much too large for an accurate estimation of the critical exponent (3 in the right- 
hand-side picture in Fig. |3 Nevertheless, we can see that the qualitative behavior, 
at least, is just such as expected from our theory. 

11 Conclusions 

Summarizing the present work we conclude the following: 

1. Critical exponents and corrections to scaling for different physical quantities 
have been discussed in framework of our [5] recently developed GFD (grouping 
of Feynman diagrams) theory (Sec. [SJ. 

2. Calculation of the two-point correlation function of 2D Ising model at the 
critical point has been made numerically by exact transfer matrix algorithms 
(Sees. |21 and 0J). The results for finite lattices including up to 800 spins have 
shown the existence of a nontrivial correction to finite-size scaling with a very 
small amplitude and exponent about 1/4, as it can be expected from our GFD 
theory. 

3. Accurate Monte Carlo simulations of the magnetization of 3D Ising model have 
been performed by Wolff's algorithm in the range of the reduced temperatures 
t > 0.000086 and system sizes L < 410 to evaluate the effective critical ex- 
ponent f3 e ff(t) based on the finite-size scaling. Estimates extracted from the 
data relatively far away from the critical point, within t > 0.0005, well confirm 
the value P ~ 0.3258 of the perturbative RG theory 47_. However, the effec- 
tive exponent tends to increase above this value when approaching T c . A self 
consistent extrapolation does not reveal a contradiction with the prediction 
P = 3/8 of the GFD theory [Hj, although there is still a large gap between 
the simulated and extrapolated values. The convergence of P e jf to the value 
of GFD theory P = 11/28 has been observed experimentally in Ni [1H], where 
the asymptotic value 0.390(3) has been found. 

4. An iterative method has been proposed (Sec. 03) which allows a direct simula- 
tion of the maximal values of the specific heat, depending on the system size L. 
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The simulated data for 3D Ising model within 6 < L < 128 apparently show a 
better agreement with the logarithmic critical singularity of the specific heat 
predicted in j)6J (and consistent with our result a = 0) than with the specific 
power-like singularity proposed by the perturbative RG theory |47j . 
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